The following tutorial is a step-by-step demo of the methods presented in our paper [1]. You may access the code and datasets here to experiment with it by yourselves.
Please make sure the following libraries are installed on your computer, as well as OpenBUGS –
library(R2OpenBUGS)
library(dplyr)
library(doParallel)
library(doSNOW)
library(tictoc)
library(tibble)
library(multinma)
library(tidyr)
library(ggplot2)
library(broom)
library(broom.mixed)
library(reshape2)
library(lemon)
library(kableExtra)
library(purrr)
library(plotly)
library(shiny)
Setting the working directory and sourcing the custom functions –
#Setting working directory to current one
setwd(dirname(rstudioapi::getActiveDocumentContext()$path))
#Sourcing functions
source('NMI_Functions.R')
The following example is one instance of a random dataset – one of thousands that were used in the simulation study in our paper. Here, we consider an evidence network consisting of four treatments, denoted \(A\), \(B\), \(C\) and \(D\), across seven studies: six – three \(A\)–\(B\) and three \(A\)–\(C\) trials – reported at the aggregate-level, and one \(A\)–\(D\) trial reported at the individual patient-level. The data were generated using the following non-shared effect modification logistic model - \[ \log\frac{p_i}{1-p_i} = \varepsilon_{i} + \begin{cases} \begin{array}{ll} -1.39 + 0.69T_i + x_{1i}T_i, &\text{$A$-$B$ trials},\\[.5em] -1.39 + T_i + 1.61x_{1i}T_i + x_{2i}T_i,&\text{$A$-$C$ trials},\\[.5em] -1.39 + 1.5T_i - 1.2x_{1i}T_i - x_{2i}T_i,&\text{$A$-$D$ trial}, \end{array} \end{cases} \] where \(x_1\) and \(x_2\) denote two binary effect modifiers, \(\varepsilon_{i}\sim\mathcal{N}\left(0,0.2^2\right)\) and \[T_i = \begin{cases} \begin{array}{ll} 0,& \text{the $i^{\mathrm{th}}$ patient received treatment $A$},\\[1em] 1,& \text{otherwise}. \end{array} \end{cases}\]
The goal of this tutorial is to demonstrate the different data requirements for the various indirect treatment comparison (ITC) methods, fit them and compare the results. In the example that will follow, the different aggregate-level datasets (AgD) were all derived from the same simulated data, and are therefore comparable. Next, we will read and display each dataset.
Starting with the IPD study –
#reading the data
IPD = read.csv('Example_IPD.csv') #IPD (for all methods)
head(IPD, 10)
## x1 x2 Tr Y Study TrtClass
## 1 1 0 D 0 7 Trt
## 2 0 1 D 1 7 Trt
## 3 1 1 D 0 7 Trt
## 4 0 0 A 0 7 Ctrl
## 5 1 0 A 1 7 Ctrl
## 6 0 0 D 0 7 Trt
## 7 1 0 A 0 7 Ctrl
## 8 0 1 A 0 7 Ctrl
## 9 1 1 D 0 7 Trt
## 10 0 0 A 0 7 Ctrl
The final \(\mathtt{TrtClass}\) column, containing the different treatment classes (in this case just \(\mathtt{Trt}\) and \(\mathtt{Ctrl}\)), is essential for the running of ML-NMR later.
It may also be interesting to estimate the Pearson correlation between the effect modifiers, as it plays an important role in NMI –
#estimating the Pearson correlation between x1 and x2
(rho_hat = cor(IPD[grep('x', colnames(IPD))]))
## x1 x2
## x1 1.0000000 0.2102674
## x2 0.2102674 1.0000000
This is the typical AgD format for standard network meta-analysis and network meta-regression [2] [3] [4] –
(NMR_AgD = read.csv('Example_AgD_NMR_NMA.csv')) #AgD for NMR/NMA
## Study Trt1 Trt2 n x1 x2 TE se
## 1 1 A B 600 0.628 0.457 1.387 0.188
## 2 2 A B 600 0.843 0.618 1.475 0.183
## 3 3 A B 600 0.787 0.592 1.612 0.187
## 4 4 A C 600 0.733 0.503 2.757 0.204
## 5 5 A C 600 0.705 0.468 2.526 0.198
## 6 6 A C 600 0.600 0.405 2.134 0.191
## 7 7 A D 600 0.533 0.303 0.667 0.187
The standard NMA will ignore the \(\mathtt{x_1}\) and \(\mathtt{x_2}\) columns, while NMR will use them as regressors.
Now for the NMI aggregate-level data in its particular format –
(NMI_AgD = read.csv('Example_AgD_NMI.csv')) #AgD for NMI
## Study Trt1 Trt2 n x1 x2 TE se
## 1 1 A B 600 0.628 0.457 1.387 0.188
## 2 1 A B 1 1.829 0.242
## 3 1 A B 0 0.492 0.319
## 4 1 A B 1 1.363 0.267
## 5 1 A B 0 1.464 0.270
## 6 2 A B 600 0.843 0.618 1.475 0.183
## 7 2 A B 1 1.577 0.198
## 8 2 A B 0 0.946 0.525
## 9 2 A B 1 1.517 0.232
## 10 2 A B 0 1.412 0.298
## 11 3 A B 600 0.787 0.592 1.612 0.187
## 12 3 A B 1 1.858 0.215
## 13 3 A B 0 0.643 0.398
## 14 3 A B 1 1.662 0.241
## 15 3 A B 0 1.542 0.296
## 16 4 A C 600 0.733 0.503 2.757 0.204
## 17 4 A C 1 3.415 0.270
## 18 4 A C 0 1.879 0.392
## 19 4 A C 1 3.135 0.308
## 20 4 A C 0 2.515 0.286
## 21 5 A C 600 0.705 0.468 2.526 0.198
## 22 5 A C 1 3.412 0.272
## 23 5 A C 0 1.038 0.329
## 24 5 A C 1 3.264 0.325
## 25 5 A C 0 2.033 0.258
## 26 6 A C 600 0.6 0.405 2.134 0.191
## 27 6 A C 1 3.139 0.279
## 28 6 A C 0 1.056 0.287
## 29 6 A C 1 3.768 0.380
## 30 6 A C 0 1.307 0.233
Note that every five rows in the above table refer to a single study, with the first row representing the study-level outcomes and the following four ar subgroup analyses. Later on, we will impute this table.
Finally, the multilevel NMR (ML-NMR) [5] aggregate-level data –
(ML_NMR_AgD = read.csv('Example_AgD_ML_NMR.csv')) #AgD for ML-NMR
## Study Tr n r X1 X2 TrtClass
## 1 1 A 318 58 0.5974843 0.4716981 Ctrl
## 2 1 B 282 133 0.6631206 0.4397163 Trt
## 3 2 A 317 65 0.8422713 0.6246057 Ctrl
## 4 2 B 283 150 0.8445230 0.6113074 Trt
## 5 3 A 292 57 0.7671233 0.5890411 Ctrl
## 6 3 B 308 169 0.8051948 0.5941558 Trt
## 7 4 A 301 56 0.7508306 0.4950166 Ctrl
## 8 4 C 299 234 0.7157191 0.5117057 Trt
## 9 5 A 280 58 0.6821429 0.4607143 Ctrl
## 10 5 C 320 245 0.7250000 0.4750000 Trt
## 11 6 A 293 57 0.6382253 0.3924915 Ctrl
## 12 6 C 307 206 0.5635179 0.4169381 Trt
Unlike the previous AgDs, this one is required to be provided as raw counts (for a dichotomous outcome) at the arm level, and not as (pre-calculated) relative treatment effects and standard errors.
The multinma package [6] makes plotting the network very convenient –
ML_NMR_data = list(IPD = IPD, AgD = ML_NMR_AgD)
net = multinma::combine_network(
set_ipd(ML_NMR_data$IPD %>%
mutate(TrC = Tr),
study = Study,
trt = Tr,
trt_class = TrtClass,
r = Y),
set_agd_arm(ML_NMR_data$AgD %>%
mutate(TrC = Tr),
study = Study,
trt = Tr,
trt_class = TrtClass,
r = r,
n = n),
trt_ref = "A")
par(mar = c(0,0,0,0))
plot(net, weight_edges = TRUE, weight_nodes = TRUE)
Now that all datasets have been read and presented, we may proceed to fit the different models.
First, there is a moderately-long list of inputs that need be supplied by the user for the program to run.
Starting with MCMC-related inputs –
N_chains = 3 #number of MCMC chains
N_iter = 1500 #number of MCMC iterations (in total)
burnin = 510 #number of MCMC samples to discard
n_int = 500 #number of MC integration points for ML-NMR
These are the column names the NMI functions require –
#AgD and IPD effect modifier column names
AgD_EM_cols = IPD_EM_cols = c('x1', 'x2')
Study_col = 'Study' #study column name
Trt_cols = c('Trt1', 'Trt2') #AgD treatment column names
TE_col = 'TE' #AgD treatment effect estimate column name
SE_col = 'se' #AgD standard error column name
IPD_Trt_col = 'Tr' #IPD treatment column name
outcome_col = 'Y' #IPD outcome column name
A vector of sample sizes for the AgD studies (by order of appearance) –
samp_sizes = rep(600, 6) #sample sizes for AgD studies
These are used in the calculation of the Best Linear Unbiased Predictor (BLUP) imputation, in the first step of NMI.
And, finally, the effect modifier levels at which we wish to perform ITC –
#value of effect modifier x1 to be used in ITC
x1 = .675
#value of effect modifier x2 to be used in ITC
x2 = .475
x_vect = c(x1, x2)
NMI uses the imputed AgD to solve two linear systems, in order to interpolate the treatment effect estimates and standard errors at \(\mathtt{\text{x_vect}}\) for all studies.
This single command performs the first two NMI steps (imputation and interpolation) –
#Imputing the AgD
NMI_object = NMI_interpolation(IPD, NMI_AgD, x_vect, AgD_EM_cols, IPD_EM_cols,
Study_col, samp_sizes, Trt_cols, TE_col, SE_col,
IPD_Trt_col)
This is the result of the first NMI step. Here, the last five rows comprise the IPD analyses.
#Have a look at the imputed dataset
NMI_object$Imputed
## Study Trt1 Trt2 x1 x2 TE se
## 1 1 A B 0.628 0.457 1.387 0.188
## 2 1 A B 1.000 0.538 1.829 0.242
## 3 1 A B 0.000 0.321 0.492 0.319
## 4 1 A B 0.739 1.000 1.363 0.267
## 5 1 A B 0.535 0.000 1.464 0.270
## 6 2 A B 0.843 0.618 1.475 0.183
## 7 2 A B 1.000 0.662 1.577 0.198
## 8 2 A B 0.000 0.381 0.946 0.525
## 9 2 A B 0.903 1.000 1.517 0.232
## 10 2 A B 0.746 0.000 1.412 0.298
## 11 3 A B 0.787 0.592 1.612 0.187
## 12 3 A B 1.000 0.646 1.858 0.215
## 13 3 A B 0.000 0.393 0.643 0.398
## 14 3 A B 0.858 1.000 1.662 0.241
## 15 3 A B 0.683 0.000 1.542 0.296
## 16 4 A C 0.733 0.503 2.757 0.204
## 17 4 A C 1.000 0.566 3.415 0.270
## 18 4 A C 0.000 0.329 1.879 0.392
## 19 4 A C 0.825 1.000 3.135 0.308
## 20 4 A C 0.639 0.000 2.515 0.286
## 21 5 A C 0.705 0.468 2.526 0.198
## 22 5 A C 1.000 0.536 3.412 0.272
## 23 5 A C 0.000 0.306 1.038 0.329
## 24 5 A C 0.807 1.000 3.264 0.325
## 25 5 A C 0.615 0.000 2.033 0.258
## 26 6 A C 0.600 0.405 2.134 0.191
## 27 6 A C 1.000 0.489 3.139 0.279
## 28 6 A C 0.000 0.279 1.056 0.287
## 29 6 A C 0.725 1.000 3.768 0.380
## 30 6 A C 0.515 0.000 1.307 0.233
## 31 7 A D 0.533 0.303 0.667 0.187
## 32 7 A D 1.000 0.394 0.008 0.278
## 33 7 A D 0.000 0.200 1.305 0.266
## 34 7 A D 0.692 1.000 -0.286 0.374
## 35 7 A D 0.464 0.000 0.997 0.223
Compare the imputed NMI AgD with the original table in the data section.
Here are the interpolated TEs and SEs, all at \(x_1 = 0.675\) and \(x_2 = 0.475\) –
#The data submitted to NMA for ITC
NMI_object$Final
## Study Trt1 Trt2 x1 x2 TE se
## 1 1 A B 0.675 0.475 1.439 0.189
## 2 2 A B 0.675 0.475 1.370 0.216
## 3 3 A B 0.675 0.475 1.480 0.196
## 4 4 A C 0.675 0.475 2.786 0.200
## 5 5 A C 0.675 0.475 2.571 0.195
## 6 6 A C 0.675 0.475 2.540 0.209
## 7 7 A D 0.675 0.475 0.289 0.221
We can visually inspect the goodness of the interpolation by plotting the originally observed treatment effect estimates against their interpolated counterparts (at the same \(x_1\) and \(x_2\) values) –
NMI_diagnostic_plotly(NMI_object)
In case one wonders about the goodness of the standard error interpolation: that system is an underdetermined one, and as such, interpolates the original “in-sample” observations without any error (meaning all point in an equivalent plot will align perfectly along the \(y=x\) line).
To conclude NMI, the final table of the previous section is submitted to standard.
#Model fitting
NMI_sim = NMA_run(NMI_object$Final, N_chains, N_iter, burnin)
#Summarising results
NMI_summary = NMA_NMI_summary(NMI_sim)
#TEs and CrIs
NMI_results = result_table(NMI_summary)
The results are stored and will eventually be displayed alongside those of the remaining methods.
#*#Model fitting
NMA_sim = NMA_run(NMR_AgD, N_chains, N_iter, burnin)
#Summarising results
NMA_summary = NMA_NMI_summary(NMA_sim)
#TEs and CrIs
NMA_results = result_table(NMA_summary)
#*#Model fitting
NMR_sim = NMA_Meta_Reg_run_2D(NMR_AgD, N_chains, N_iter, burnin)
#Summarising results
NMR_summary = NMA_Metareg_summary_2D(NMR_sim, x_vect)
#TEs and CrIs
NMR_results = result_table(NMR_summary)
#Model fitting
ML_NMR_sim = ML_NMR_Run_2D(ML_NMR_data, N_iter, N_chains, burnin, n_int)
#summarising results
ML_NMR_summ = ML_NMR_summary_2D(n_trts = 4, ML_NMR_sim, x_vect)
#TEs and CrIs
ML_NMR_results = result_table(ML_NMR_summ)
Let us now compare the results of the various methods, relative to true parameters used in the data generation.
We start by calculating the true (relative) treatment effects from the logistic model parameters -
#These are the logistic regression coefficients that were used to generate the
#data used in this demonstration
beta_AB = c(-1.39, 0, 0, 0.69, 1.00, 0.00)
beta_AC = c(-1.39, 0, 0, 1.00, 1.61, 1.00)
beta_AD = c(-1.39, 0, 0, 1.50, -1.20, -1.00)
#converting the logistic regression coefficients into treatment effects
d = rbind(beta_AB, beta_AC, beta_AD)[,4:6]
trt_effs = d%*%c(1, x1, x2)
(trt_effs = c(trt_effs, trt_effs[c(2,3,3)] - trt_effs[c(1,1,2)]))
## [1] 1.36500 2.56175 0.21500 1.19675 -1.15000 -2.34675
We can now put together everything we have done.
display_result_table(NMA_results, NMR_results, ML_NMR_results, NMI_results, trt_effs)
| NMA | NMR | ML-NMR | NMI | True Trt. Eff. | |
|---|---|---|---|---|---|
| \(d_{\mathrm{AB}}\) | 1.49 [1.28 ;1.70] | 1.46 [1.17 ;1.73] | 1.54 [1.32 ;1.75] | 1.43 [1.21 ;1.65] | 1.36 |
| \(d_{\mathrm{AC}}\) | 2.46 [2.24 ;2.68] | 2.38 [2.12 ;2.61] | 2.66 [2.40 ;2.91] | 2.63 [2.41 ;2.86] | 2.56 |
| \(d_{\mathrm{AD}}\) | 0.67 [0.30 ;1.02] | 0.61 [0.07 ;1.19] | 0.21 [-0.24 ;0.73] | 0.29 [-0.15 ;0.70] | 0.21 |
| \(d_{\mathrm{BC}}\) | 0.96 [0.66 ;1.29] | 0.92 [0.52 ;1.33] | 1.12 [0.79 ;1.43] | 1.20 [0.88 ;1.54] | 1.20 |
| \(d_{\mathrm{BD}}\) | -0.83 [-1.26 ;-0.41] | -0.86 [-1.56 ;-0.07] | -1.32 [-1.81 ;-0.77] | -1.15 [-1.64 ;-0.67] | -1.15 |
| \(d_{\mathrm{CD}}\) | -1.79 [-2.22 ;-1.39] | -1.76 [-2.30 ;-1.18] | -2.44 [-3.00 ;-1.86] | -2.35 [-2.84 ;-1.89] | -2.35 |
result_forest_plot(NMA_summary, NMR_summary, ML_NMR_summ, NMI_summary, trt_effs)
## Warning: Using `size` aesthetic for lines was deprecated in ggplot2 3.4.0.
## ℹ Please use `linewidth` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
In this example, NMI did extremely well in terms of estimation accuracy (especially in the case of the indirect comparisons, namely: \(d_{\mathrm{BC}}\), \(d_{\mathrm{BD}}\) and \(d_{\mathrm{CD}}\)). Note that while this tutorial compares various ITC methods on a single random dataset, the paper itself contains results from many thousands of simulated datasets. We strongly encourage you to give it a read before drawing any conclusions.
NMI is an attractive option when subgroup analyses for all effect modifiers are available. The extra work put toward data extraction, more than pays itself off when the data is driven by a non-shared effect modification data generating mechanism.